#include <assert.h>
#include <stdio.h>

#include "../src/mstep.c"

static void
fill (double*           x,
      double          (*func) (void),
      unsigned int      n)
{
  for (unsigned int i = 0; i < n; ++i)
    { 
      x[i] = func ();
    }
}

static void
fill_r (Rating*         x,
        unsigned int    n_workers,
        unsigned int    n_labels,
        unsigned int    n_images,
        unsigned int    n)
{
  for (unsigned int i = 0; i < n; ++i)
    { 
      x[i].worker = (lrand48 () >> 16) % n_workers;
      x[i].image = (lrand48 () >> 16) % n_images;
      x[i].rating = (lrand48 () >> 16) % n_labels;
    }
}

static void
test_to_from_identity_once (void)
{ 
  /* (kappa, logalpha, tau, gamma) = 
     from_x (kappa, logalpha, tau, gamma, to_x (kappa, logalpha, tau, gamma)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  double kappa;
  double logalpha[n_images];
  double tau[n_workers * (n_labels - 1)];
  double gamma[n_labels - 1];
  double rho;
  double loglambda[n_workers];

  double newkappa;
  double newlogalpha[n_images];
  double newtau[n_workers * (n_labels - 1)];
  double newgamma[n_labels - 1];
  double newrho;
  double newloglambda[n_workers];

  kappa = drand48 ();
  fill (logalpha, drand48, n_images);
  fill (tau, drand48, n_workers * (n_labels - 1));
  fill (gamma, drand48, n_labels - 1);
  rho = drand48 ();
  fill (loglambda, drand48, n_workers);

  newkappa = drand48 ();
  fill (newlogalpha, drand48, n_images);
  fill (newtau, drand48, n_workers * (n_labels - 1));
  fill (newgamma, drand48, n_labels - 1);
  newrho = drand48 ();
  fill (newloglambda, drand48, n_workers);

  size_t n_dims = 
    1 + n_images + n_workers * (n_labels - 1) + (n_labels - 1) + 1 + n_workers;
  gsl_vector* x = gsl_vector_alloc (n_dims);

  to_x (&kappa, logalpha, n_workers, n_labels, tau, n_images, gamma, &rho, loglambda, x);
  from_x (&newkappa, newlogalpha, n_workers, n_labels, newtau, n_images, newgamma, &newrho, newloglambda, x);
          
  gsl_vector_free (x);

  assert (kappa == newkappa);
  assert (memcmp (logalpha, newlogalpha, n_images * sizeof (logalpha[0])) == 0);
  assert (memcmp (tau,
                  newtau,
                  n_workers * (n_labels - 1) * sizeof (tau[0])) == 0);
  assert (memcmp (gamma, newgamma, (n_labels - 1) * sizeof (gamma[0])) == 0);
  assert (rho == newrho);
  assert (memcmp (loglambda, newloglambda, n_workers * sizeof (loglambda[0])) == 0);
}

static void
test_to_from_inverse (void)
{
  for (unsigned int i = 0; i < 100; ++i)
    { 
      test_to_from_identity_once ();
    }
}

static void
test_qfunc_f_once (void)
{
  /* qfunc (kappa, logalpha, tau, gamma) = 
     -qfunc_f (to_x (kappa, logalpha, tau, gamma)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_ratings = ((lrand48 () >> 16) % 5) + 1;

  double kappa;
  double logalpha[n_images];
  double tau[n_workers * (n_labels - 1)];
  double gamma[n_labels - 1];
  double rho;
  double loglambda[n_workers];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  Rating ratings[n_ratings];
  NormalDistribution priorkappa = normal_distribution (1.0 + drand48 (),
                                                       1.0 + drand48 ());
  ParameterizedNormalDistribution priorlogalpha = 
    parameterized_normal_distribution (1.0 + drand48 ());
  ParameterizedNormalDistribution priortau = 
    parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution normalpriorgamma[n_labels - 1];
  Distribution* priorgamma[n_labels - 1];
  NormalDistribution priorrho = normal_distribution (1.0 + drand48 (),
                                                     1.0 + drand48 ());
  ParameterizedNormalDistribution priorloglambda =
    parameterized_normal_distribution (1.0 + drand48 ());

  for (unsigned int i = 0; i + 1 < n_labels; ++i)
    {
      normalpriorgamma[i] = normal_distribution (1.0 + drand48 (),
                                                 1.0 + drand48 ());
      priorgamma[i] = &normalpriorgamma[i].base;
    }

  kappa = drand48 ();
  fill (logalpha, drand48, n_images);
  fill (tau, drand48, n_workers * (n_labels - 1));
  fill (gamma, drand48, n_labels - 1);
  rho = drand48 ();
  fill (loglambda, drand48, n_workers);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);
  fill_r (ratings, n_workers, n_labels, n_images, n_ratings);

  QfuncParams qfunc_params = { .n_workers = n_workers,
                               .n_labels = n_labels,
                               .n_images = n_images,
                               .pz = pz,
                               .priorz = priorz,
                               .ratings = ratings,
                               .n_ratings = n_ratings,
                               .priorkappa = &priorkappa.base,
                               .priorlogalpha = &priorlogalpha.base,
                               .priortau = &priortau.base,
                               .priorgamma = priorgamma,
                               .priorrho = &priorrho.base,
                               .priorloglambda = &priorloglambda.base };

  size_t n_dims = 
    1 + n_images + n_workers * (n_labels - 1) + (n_labels - 1) + 1 + n_workers;
  gsl_vector* x = gsl_vector_alloc (n_dims);

  to_x (&kappa, logalpha, n_workers, n_labels, tau, n_images, gamma, &rho, loglambda, x);

  double q = qfunc (kappa, logalpha, n_workers, n_labels,
                    tau, n_images, gamma, rho, loglambda, pz, priorz,
                    ratings, n_ratings, &priorkappa.base,
                    &priorlogalpha.base, &priortau.base,
                    priorgamma, &priorrho.base, &priorloglambda.base);

  double qf = qfunc_f (x, &qfunc_params);

  assert (fabs (q + qf) <= 1e-8 * (1 + fabs (q) + fabs (qf)) ||
          (fprintf (stderr, "%g ?= %g\n", q, qf), 0));

  gsl_vector_free (x);
}

static void
test_qfunc_f (void)
{ 
  for (unsigned int i = 0; i < 100; ++i)
    { 
      test_qfunc_f_once ();
    }
}

static void
test_qfunc_df_once (void)
{
  /* dqfunc (kappa, logalpha, tau, gamma) = 
     -qfunc_df (to_x (kappa, logalpha, tau, gamma)) */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_ratings = ((lrand48 () >> 16) % 5) + 1;

  double kappa;
  double logalpha[n_images];
  double tau[n_workers * (n_labels - 1)];
  double gamma[n_labels - 1];
  double rho;
  double loglambda[n_workers];
  double dkappa;
  double dlogalpha[n_images];
  double dtau[n_workers * (n_labels - 1)];
  double dgamma[n_labels - 1];
  double drho;
  double dloglambda[n_workers];
  double newdkappa;
  double newdlogalpha[n_images];
  double newdtau[n_workers * (n_labels - 1)];
  double newdgamma[n_labels - 1];
  double newdrho;
  double newdloglambda[n_workers];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  Rating ratings[n_ratings];
  NormalDistribution priorkappa = normal_distribution (1.0 + drand48 (),
                                                       1.0 + drand48 ());
  ParameterizedNormalDistribution priorlogalpha = 
    parameterized_normal_distribution (1.0 + drand48 ());
  ParameterizedNormalDistribution priortau = 
    parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution normalpriorgamma[n_labels - 1];
  Distribution* priorgamma[n_labels - 1];
  NormalDistribution priorrho = normal_distribution (1.0 + drand48 (),
                                                     1.0 + drand48 ());
  ParameterizedNormalDistribution priorloglambda =
    parameterized_normal_distribution (1.0 + drand48 ());
    
  for (unsigned int i = 0; i + 1 < n_labels; ++i)
    {
      normalpriorgamma[i] = normal_distribution (1.0 + drand48 (),
                                                 1.0 + drand48 ());
      priorgamma[i] = &normalpriorgamma[i].base;
    }

  kappa = drand48 ();
  fill (logalpha, drand48, n_images);
  fill (tau, drand48, n_workers * (n_labels - 1));
  fill (gamma, drand48, n_labels - 1);
  rho = drand48 ();
  fill (loglambda, drand48, n_workers);
  dkappa = drand48 ();
  fill (dlogalpha, drand48, n_images);
  fill (dtau, drand48, n_workers * (n_labels - 1));
  fill (dgamma, drand48, n_labels - 1);
  drho = drand48 ();
  fill (dloglambda, drand48, n_workers);
  newdkappa = drand48 ();
  fill (newdlogalpha, drand48, n_images);
  fill (newdtau, drand48, n_workers * (n_labels - 1));
  fill (newdgamma, drand48, n_labels - 1);
  newdrho = drand48 ();
  fill (newdloglambda, drand48, n_workers);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);
  fill_r (ratings, n_workers, n_labels, n_images, n_ratings);

  QfuncParams qfunc_params = { .n_workers = n_workers,
                               .n_labels = n_labels,
                               .n_images = n_images,
                               .pz = pz,
                               .priorz = priorz,
                               .ratings = ratings,
                               .n_ratings = n_ratings,
                               .priorkappa = &priorkappa.base,
                               .priorlogalpha = &priorlogalpha.base,
                               .priortau = &priortau.base,
                               .priorgamma = priorgamma,
                               .priorrho = &priorrho.base,
                               .priorloglambda = &priorloglambda.base };

  size_t n_dims = 
    1 + n_images + n_workers * (n_labels - 1) + (n_labels - 1) + 1 + n_workers;
  gsl_vector* x = gsl_vector_alloc (n_dims);
  gsl_vector* g = gsl_vector_alloc (n_dims);

  to_x (&kappa, logalpha, n_workers, n_labels, tau, n_images, gamma, &rho, loglambda, x);

  dqfunc (&dkappa, dlogalpha, dtau, dgamma, &drho, dloglambda,
          kappa, logalpha, n_workers, n_labels,
          tau, n_images, gamma, rho, loglambda, 
          pz, ratings, n_ratings, &priorkappa.base,
          &priorlogalpha.base, &priortau.base,
          priorgamma, &priorrho.base, &priorloglambda.base);

  qfunc_df (x, &qfunc_params, g);

  from_x (&newdkappa, newdlogalpha, n_workers, n_labels, newdtau, n_images, newdgamma, &newdrho, newdloglambda, g);

  assert (fabs (dkappa + newdkappa) <= 
          1e-8 * (1 + fabs (dkappa) + fabs (newdkappa)) ||
          (fprintf (stderr, "%g ?= %g\n", dkappa, newdkappa), 0));

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (fabs (dlogalpha[i] + newdlogalpha[i]) <=
              1e-8 * (1 + fabs (dlogalpha[i]) + fabs (newdlogalpha[i])) ||
              (fprintf (stderr, "%g ?= %g\n", dlogalpha[i], newdlogalpha[i]), 0));
    }

  for (unsigned int i = 0; i < n_workers * (n_labels - 1); ++i)
    {
      assert (fabs (dtau[i] + newdtau[i]) <=
              1e-8 * (1 + fabs (dtau[i]) + fabs (newdtau[i])) ||
              (fprintf (stderr, "%g ?= %g\n", dtau[i], newdtau[i]), 0));
    }

  for (unsigned int i = 0; i < (n_labels - 1); ++i)
    {
      assert (fabs (dgamma[i] + newdgamma[i]) <=
              1e-8 * (1 + fabs (dgamma[i]) + fabs (newdgamma[i])) ||
              (fprintf (stderr, "%g ?= %g\n", dgamma[i], newdgamma[i]), 0));
    }

  assert (fabs (drho + newdrho) <= 
          1e-8 * (1 + fabs (drho) + fabs (newdrho)) ||
          (fprintf (stderr, "%g ?= %g\n", drho, newdrho), 0));

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      assert (fabs (dloglambda[i] + newdloglambda[i]) <=
              1e-8 * (1 + fabs (dloglambda[i]) + fabs (newdloglambda[i])) ||
              (fprintf (stderr, "%g ?= %g\n", dloglambda[i], newdloglambda[i]), 0));
    }

  gsl_vector_free (g);
  gsl_vector_free (x);
}

static void
test_qfunc_df (void)
{ 
  for (unsigned int i = 0; i < 100; ++i)
    { 
      test_qfunc_df_once ();
    }
}

static void
test_mstep_no_data_once (void)
{
  /* without data mstep should converge to prior maximum */

  unsigned int n_workers = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_labels = ((lrand48 () >> 16) % 5) + 1;
  unsigned int n_images = ((lrand48 () >> 16) % 5) + 1;

  double kappa;
  double logalpha[n_images];
  double tau[n_workers * (n_labels - 1)];
  double gamma[n_labels - 1];
  double rho;
  double loglambda[n_workers];
  double pz[n_images * n_labels];
  double priorz[n_labels];
  NormalDistribution priorkappa = normal_distribution (1.0 + drand48 (),
                                                       1.0 + drand48 ());
  ParameterizedNormalDistribution priorlogalpha = 
    parameterized_normal_distribution (1.0 + drand48 ());
  ParameterizedNormalDistribution priortau = 
    parameterized_normal_distribution (1.0 + drand48 ());
  NormalDistribution normalpriorgamma[n_labels - 1];
  Distribution* priorgamma[n_labels - 1];
  NormalDistribution priorrho = normal_distribution (1.0 + drand48 (),
                                                     1.0 + drand48 ());
  ParameterizedNormalDistribution priorloglambda =
    parameterized_normal_distribution (1.0 + drand48 ());

  for (unsigned int i = 0; i + 1 < n_labels; ++i)
    {
      normalpriorgamma[i] = normal_distribution (1.0 + drand48 (),
                                                 1.0 + drand48 ());
      priorgamma[i] = &normalpriorgamma[i].base;
    }

  kappa = drand48 ();
  fill (logalpha, drand48, n_images);
  fill (tau, drand48, n_workers * (n_labels - 1));
  fill (gamma, drand48, n_labels - 1);
  rho = drand48 ();
  fill (loglambda, drand48, n_workers);
  fill (pz, drand48, n_images * n_labels);
  fill (priorz, drand48, n_labels);

  mstep (&kappa, logalpha, n_workers, n_labels,
         tau, n_images, gamma, &rho, loglambda,
         pz, priorz, NULL, 0, 
         &priorkappa.base, &priorlogalpha.base,
         &priortau.base, priorgamma, &priorrho.base, &priorloglambda.base,
         100, 1e-8);

  assert (   fabs (kappa - priorkappa.mean) 
          <= 1e-4 * (1 + fabs (kappa) + fabs (priorkappa.mean))
          || (fprintf (stderr, "%g ?= %g\n", kappa, priorkappa.mean), 0));

  for (unsigned int i = 0; i < n_images; ++i)
    {
      assert (   fabs (logalpha[i] - priorkappa.mean)
              <= 1e-4 * (1 + fabs (logalpha[i]) + fabs (priorkappa.mean))
              || (fprintf (stderr, "%u: %g ?= %g\n", i, logalpha[i], priorkappa.mean), 0));
    }

  for (unsigned int i = 0; i + 1 < n_labels; ++i)
    {
      assert (   fabs (gamma[i] - normalpriorgamma[i].mean)
              <= 1e-4 * (1 + fabs (gamma[i]) + fabs (normalpriorgamma[i].mean))
              || (fprintf (stderr, "%u: %g ?= %g\n", i, gamma[i], normalpriorgamma[i].mean), 0));
    }

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      for (unsigned int j = 0; j + 1 < n_labels; ++j)
        {
          assert (   fabs (tau[i * (n_labels - 1) + j] - normalpriorgamma[j].mean)
                  <= 1e-4 * (1 + fabs (tau[i * (n_labels - 1) + j]) + fabs (normalpriorgamma[j].mean))
                  || (fprintf (stderr, "%u: %g ?= %g\n", i, tau[i * (n_labels - 1) + j], normalpriorgamma[j].mean), 0));
        }
    }

  assert (   fabs (rho - priorrho.mean) 
          <= 1e-4 * (1 + fabs (rho) + fabs (priorrho.mean))
          || (fprintf (stderr, "%g ?= %g\n", rho, priorrho.mean), 0));

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      assert (   fabs (loglambda[i] - priorrho.mean)
              <= 1e-4 * (1 + fabs (loglambda[i]) + fabs (priorrho.mean))
              || (fprintf (stderr, "%u: %g ?= %g\n", i, loglambda[i], priorrho.mean), 0));
    }
}

static void
test_mstep_no_data (void)
{
  for (unsigned int i = 0; i < 20; ++i)
    {
      test_mstep_no_data_once ();
    }
}

typedef struct _MstepTestDatum MstepTestDatum;
struct _MstepTestDatum 
{
  struct
    {
      double kappa;
      double logalpha[3];
      double tau[4*(5-1)];
      double gamma[(5-1)];
      double rho;
      double loglambda[4];
      double pz[3*5];
      double priorz[5];
      Rating ratings[8];
      double priorlogalphastddev;
      double priortaustddev;
      double priorgamma1mean;
      double priorgamma1stddev;
      double priorgamma2mean;
      double priorgamma2stddev;
      double priorgamma3mean;
      double priorgamma3stddev;
      double priorgamma4mean;
      double priorgamma4stddev;
      double priorkappamean;
      double priorkappastddev;
      double priorrhomean;
      double priorrhostddev;
      double priorloglambdastddev;
    }                           input;

  struct
    {
      double q;
      struct
        {
          double kappa;
          double logalpha[3];
          double tau[4*(5-1)];
          double gamma[(5-1)];
          double rho;
          double loglambda[4];
        }    v;
    }                           output;
};

static void
test_mstep (void)
{
  MstepTestDatum data[] = 
    {{{0.469458, {0.415694, 0.00293394, 0.33327}, {0.0174874, 0.900134,  0.146335, 0.00937511, 0.194021, 0.924243, 0.769001, 0.521991,  0.354686, 0.383854, 0.484947, 0.86034, 0.897873, 0.319177,  0.350202, 0.240425}, {0.411082, 0.267254, 0.742008, 0.206705},  0.941625, {0.85156, 0.739074, 0.873435, 0.924138}, {0.951425,  0.592739, 0.86406, 0.730117, 0.0271821, 0.823737, 0.342068,  0.375431, 0.643328, 0.33879, 0.481729, 0.477558, 0.324151,  0.988589, 0.241304}, {0.0664752, 0.0568969, 0.246581, 0.0345992,  0.12485}, {{0, 2, 4}, {1, 1, 1}, {3, 1, 3}, {2, 0, 1}, {2, 1,  2}, {1, 1, 2}, {0, 1, 0}, {2, 1, 4}}, 0.455337, 0.757507,  0.161165, 0.450713, 0.253912, 1.16477, 0.297105, 0.720596, 0.22673,  0.34103, 0.955037, 0.345165, 0.583402, 1.00224,  0.723308}, {-75.039476107536471083, {0.68027549356845052743, {0.48577544532116316960, 0.19960770932314396609,  0.87728840599525799996}, {0.30283236950707197288,  0.73033104921040698625, 0.44792384128258457789,  0.051703326780473268414, 0.088979563726891206426,  0.73080708839130528959, 0.67137017474126224794,  0.33964025624417347143, 0.035759507083079114466,  0.77409802740732308215, 0.66306752136104695397,  0.24611587575580215992, 0.10783346841702818497,  0.57872866644898512972, 0.34820311202288620740,  0.33435649854822017359}, {0.14515603590240888507,  0.66049919398195168040, 0.48165579402113008202,  0.23399404605956341801}, -1.7331690220513349288, {-2.0317378458536190646, -2.2202348285077705260, -2.1221222756535828916, -1.7651411752081027452}}}}, {{0.617607, {0.259252, 0.763651, 0.232004}, {0.551132, 0.202355, 0.517071,  0.197405, 0.426282, 0.997017, 0.00956392, 0.0362408, 0.225569,  0.743105, 0.0947962, 0.739136, 0.754973, 0.516375, 0.00376595,  0.784099}, {0.659808, 0.932973, 0.251526, 0.310791},  0.0422012, {0.673722, 0.487875, 0.0787865, 0.491069}, {0.471367,  0.970804, 0.881381, 0.0647876, 0.47435, 0.96124, 0.84514,  0.839219, 0.731244, 0.866444, 0.106005, 0.0842454, 0.214869,  0.862678, 0.321905}, {0.424437, 0.281896, 0.611152, 0.0111145,  0.382236}, {{1, 1, 1}, {0, 2, 0}, {2, 0, 4}, {2, 2, 1}, {0, 0,  3}, {0, 2, 3}, {3, 0, 0}, {3, 2, 4}}, 0.858174, 0.373277,  0.932328, 1.14117, 0.136808, 0.402473, 0.0509468, 1.07638,  0.662458, 0.441233, 0.205806, 1.23716, 0.931214, 0.574789,  0.349802}, {-52.697572161500072472, {-2.2699327671409699373, {-2.3584558277094818637, -2.8712273954595612360, -2.7713686620716040031}, {1.2121346706236245581,  0.18742340715909663179, -0.0018086796758328864528,  0.69389288793295426707, 1.1983438061754558999,  0.22529897398773694384, 0.044806187847708405838,  0.67867460183165053996, 1.1945566599568997333,  0.19327065670145367778, 0.0071784307423506431109,  0.63283811235109057315, 1.2472841102812357201,  0.22577422591363782826, 0.040670172805077659199,  0.66797753913991322526}, {1.2057656612457658163,  0.19535209606087200995, 0.023535667294396754851,  0.66745220354602472410},  0.39680979176744136356, {0.35952191129557674108,  0.29690632783732629636, 0.40801641464736832333,  0.32487097687047375454}}}}, {{0.902915, {0.716345, 0.462111,  0.777896}, {0.478478, 0.434449, 0.850958, 0.766782, 0.096242,  0.826274, 0.727681, 0.834454, 0.205076, 0.689467, 0.575208,  0.783507, 0.378697, 0.0270088, 0.383975, 0.577701}, {0.391537,  0.0957951, 0.059186, 0.477899},  0.488622, {0.37945, 0.597075, 0.700002, 0.0101442}, {0.945002,  0.746117, 0.93322, 0.913902, 0.118727, 0.0184361, 0.0987663,  0.708827, 0.429261, 0.443228, 0.315259, 0.33013, 0.402252,  0.0592536, 0.737558}, {0.938593, 0.306457, 0.0000675239, 0.25966,  0.449971}, {{2, 2, 4}, {1, 1, 4}, {1, 0, 1}, {2, 2, 2}, {2, 2,  4}, {0, 1, 0}, {0, 2, 1}, {1, 2, 0}}, 1.17701, 0.652992,  0.559657, 0.689827, 0.982005, 0.906875, 0.626437, 0.775925,  0.863277, 0.888439, 0.527671, 1.0671, 0.434017, 0.44521,  0.462411}, {-68.847652468437792625, {-1.1942160278987720504, {-1.8349556568758088061, -1.7654376176673588726, -2.0771076856019859074}, {0.70112927461084529179,  1.5288667294991237773, 0.94384485588245785099,  0.78289509338445939423, 0.64959801764731235609,  1.5917243018617148558, 0.95565141519611964246,  0.73516311930561164056, 0.57200226307092745335,  1.2881579233306399577, 0.79795195501733567334,  0.59434878298570333563, 0.62222248348824129725,  1.3977354595000650256, 0.84706415127997282845,  0.72841943218178513788}, {0.62222248348625251864,  1.3977354595012369304, 0.84706415127993252189,  0.72841943218341761400},  0.13428927098503571296, {-0.17158326376371280839, -0.00024294607043442304339, 0.25135877709161642002,  0.13428927098501742593}}}}};

  for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
    {
      ParameterizedNormalDistribution priorlogalpha = 
        parameterized_normal_distribution (data[i].input.priorlogalphastddev);

      ParameterizedNormalDistribution priortau =
        parameterized_normal_distribution (data[i].input.priortaustddev);

      NormalDistribution priorgammanormal[4] = {
          normal_distribution (data[i].input.priorgamma1mean,
                               data[i].input.priorgamma1stddev),
          normal_distribution (data[i].input.priorgamma2mean,
                               data[i].input.priorgamma2stddev),
          normal_distribution (data[i].input.priorgamma3mean,
                               data[i].input.priorgamma3stddev),
          normal_distribution (data[i].input.priorgamma4mean,
                               data[i].input.priorgamma4stddev)
      };

      Distribution* priorgamma[4] = {
          &priorgammanormal[0].base,
          &priorgammanormal[1].base,
          &priorgammanormal[2].base,
          &priorgammanormal[3].base
      };

      NormalDistribution priorkappa = 
        normal_distribution (data[i].input.priorkappamean,
                             data[i].input.priorkappastddev);

      NormalDistribution priorrho = 
        normal_distribution (data[i].input.priorrhomean,
                             data[i].input.priorrhostddev);

      ParameterizedNormalDistribution priorloglambda =
        parameterized_normal_distribution (data[i].input.priorloglambdastddev);

      double kappa = data[i].input.kappa;
      double logalpha[3];
      double tau[4*(5-1)];
      double gamma[(5-1)];
      double rho = data[i].input.rho;
      double loglambda[4];

      memcpy (logalpha, data[i].input.logalpha, 3 * sizeof (double));
      memcpy (tau, data[i].input.tau, 4*(5-1) * sizeof (double));
      memcpy (gamma, data[i].input.gamma, (5-1) * sizeof (double));
      memcpy (loglambda, data[i].input.loglambda, 4 * sizeof (double));

      double q = mstep (&kappa,
                        logalpha,
                        4,
                        5,
                        tau,
                        3,
                        gamma,
                        &rho,
                        loglambda,
                        data[i].input.pz,
                        data[i].input.priorz,
                        data[i].input.ratings,
                        8,
                        &priorkappa.base,
                        &priorlogalpha.base,
                        &priortau.base,
                        priorgamma,
                        &priorrho.base,
                        &priorloglambda.base,
                        10000,
                        1e-16);

      assert (      fabs (q - data[i].output.q) 
                 <= 1e-4 * (1 + fabs (q) + fabs (data[i].output.q))
              || (fprintf (stderr, "fail: %g != %g\n", q, data[i].output.q),
                  0));

      assert (   fabs (kappa - data[i].output.v.kappa)
              <= 1e-4 * (1 + fabs (kappa) + fabs (data[i].output.v.kappa))
              || (fprintf (stderr, "%g ?= %g\n", kappa, data[i].output.v.kappa), 0));

      for (unsigned int j = 0; j < 3; ++j)
        {
          assert (   fabs (logalpha[j] - data[i].output.v.logalpha[j])
                  <= 1e-4 * (1 + fabs (logalpha[j]) + fabs (data[i].output.v.logalpha[j]))
                  || (fprintf (stderr, "%u: %g ?= %g\n", j, logalpha[j], data[i].output.v.logalpha[j]), 0));
        }

      for (unsigned int j = 0; j < 4; ++j)
        {
          for (unsigned int k = 0; k + 1 < 5; ++k)
            {
              assert (   fabs (tau[j * (5 - 1) + k] - data[i].output.v.tau[j * (5 - 1) + k])
                      <= 1e-4 * (1 + fabs (tau[j * (5 - 1) + k]) + fabs (data[i].output.v.tau[j * (5 - 1) + k]))
                      || (fprintf (stderr, "%u: %g ?= %g\n", j, tau[j * (5 - 1) + k], data[i].output.v.tau[j * (5 - 1) + k]), 0));
            }
        }

      for (unsigned int j = 0; j + 1 < 5; ++j)
        {
          assert (   fabs (gamma[j] - data[i].output.v.gamma[j])
                  <= 1e-4 * (1 + fabs (gamma[j]) + fabs (data[i].output.v.gamma[j]))
                  || (fprintf (stderr, "%u: %g ?= %g\n", j, gamma[j], data[i].output.v.gamma[j]), 0));
        }

      assert (   fabs (rho - data[i].output.v.rho)
              <= 1e-4 * (1 + fabs (rho) + fabs (data[i].output.v.rho))
              || (fprintf (stderr, "%g ?= %g\n", rho, data[i].output.v.rho), 0));

      for (unsigned int j = 0; j < 4; ++j)
        {
          assert (   fabs (loglambda[j] - data[i].output.v.loglambda[j])
                  <= 1e-4 * (1 + fabs (loglambda[j]) + fabs (data[i].output.v.loglambda[j]))
                  || (fprintf (stderr, "%u: %g ?= %g\n", j, loglambda[j], data[i].output.v.loglambda[j]), 0));
        }
    }
}

int
main (void)
{
  srand48 (69);

  test_to_from_inverse ();
  test_qfunc_f ();
  test_qfunc_df ();
  test_mstep_no_data ();
  test_mstep ();

  return 0;
}
